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Abstract. Time-dependent density matrix renormalization group method with a 
matrix product ansatz is employed for explicit computation of non-equilibrium steady 
state density operators of several integrable and non-integrable quantum spin chains, 
which are driven far from equilibrium by means of Markovian couplings to external 
baths at the two ends. It is argued that even though the time-evolution can not 
be simulated efficiently due to fast entanglement growth, the steady states in and 
out of equilibrium can be typically accurately approximated, so that chains of length 
of the order n Ki 100 are accessible. Our results are demonstrated by performing 
expHcit simulations of steady states and calculations of energy/spin densities/currents 
in several problems of heat and spin transport in quantum spin chains. Previously 
conjectured relation between quantum chaos and normal transport is re-confirmed 
with high acurracy on much larger systems. 



Submitted to: J. Stat. Mech. 
1. Introduction 

A detailed understanding of the properties of strongly interacting many-particle 
quantum systems represents one of the major challenges of theoretical physics. Even 
for one spatial dimension (ID) and local interactions very little is known about physics 
at high temperature, out of equilibrium, or for long real time evolution. 

Quite recently, several original ideas emerging from the quantum information theory 
on the issue of entanglement in many-particle systems [Ij resulted in very useful efficient 
techniques for the "classical" simulation of quantum many particle systems |2l [3] . Even 
though most versions of these methods can be re-interpreted in terms of the density 
matrix renormalization group |1] (for an excellent review see [5]) and its time-dependent 
generalizations [HI [7] they bring up a conceptual simplicity: namely, they represent 
the quantum many body states in terms of the so-called finitely correlated states [8] 
which can be written in terms of matrix products, the so-called matrix product states 
(MPS). Similarly, matrix product operator (MPO) ansatz can be used to describe density 
operators of mixed states [9l [10] . 
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Yet, most of these ingenious numerical methods turned out to be effective only 
for calculation of ground states, low temperature equilibrium states, or time evolved 
states for short time or starting from particular (e.g. few "quasi-particle" ) initial 
conditions. The calculation of asymptotic time evolution in the thermodynamic limit, 
and the properties of non- equilibrium steady states (NESS), remains notoriously difficult 
due to typical growth of bi-partite entanglement in generic non-integrable systems and 
consequent exponential growth of the dimension of the relevant many-particle Hilbert 
space [H]. Therefore, the regularity to quantum chaos transition is manifested in the 
difficulty of classical simulation of quantum many-body dynamics [ri|, see also |12j . 
Even though promising analytical techniques based on Fock spaces of operators have 
been developed recently in terms of which one can access NESS and relaxation properties 
of certain open integrable quantum many-body systems far from equilibrium [131 fT4] . 
for generic non-integrable systems one still needs to rely on more or less sophisticated 
"brute force" methods limited at present to 20-30 spins 1/2 (or qubits, fermions) in 
ID. For such small systems it is often difficult or impossible to conclude on asymptotic 
thermodynamic properties due to considerable finite size effects. 

In this paper we shall put forward a hypothesis that NESS of certain open ID 
systems can be efficiently simulated in terms of MPO ansatz. 

Our idea is based on the following arguments. It has been shown rigorously that 
ground states of several interacting many-particles systems in ID can be written in 
terms of MPS (for example, the so-called valence bond states [15]). Later, it has 
been shown [HI [17], using the arguments of conformal field theory, that ground state 
entanglement entropy of a finite block of size n in any non-critical (gaped) system in ID 
saturates with n, whereas for critical (non-gaped) systems the growth of entanglement 
entropy is at most logarithmic in n. This result immediately implies that an efficient 
MPS representation of ground states in ID should exist with the dimension of the 
parameterizing matrices which saturates or grows not faster than polynomially with the 
system size n. Similar "strong efficiency" result can be claimed for thermal states |18j . 
The defining equation for NESS of an open-quantum system can be written as "right 
ground state" of a certain non-Hermitian quantum Liouville master operator [13] , hence 
NESS can be thought of as a kind of non-normal ground state in the Liouville space. 
This means that the entanglement of NESS, treated as an element of the Hilbert space 
of operators, can be relatively weak for a a wide class of (non-integrable) problems even 
though the operator space entanglement of generic operators is much stronger (e.g. 
of the "excited" eigenstates of Liouville master operator, i.e. decay modes). If this 
intuition is correct then NESS should be simulable in terms of MPO ansatz for systems 
considerably longer than 30 spins (at present). 

In the present work we are going to test MPO method on NESS calculation of 
different ID systems, integrable and chaotic, as well as normal and ideal conductors. In 
particular, we shall focus on the problem of diffusive versus ballistic spin/heat transport 
in spin chains. There have been basically two approaches to quantum transport. The 
first one is using the linear response formalism, calculating equilibrium time correlation 
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functions. In that approach one studies purely Hamiltonian system, i.e., without any 
external baths, either directly (see e.g. |l9l[20l|2ll[22l[23l|2ll|25l[26]), orin terms of 
conformal field theory p^, or quantum Monte Carlo [28] • For exact calculations spin 
chains of sizes up to n = 28 are achievable using a micro-canonical finite temperature 
Lanczos method [2^. For a discussion of applicability of Green-Kubo formalism see |29j . 
The second approach goes via direct simulation of out-of-equilibrium system coupled to 
the baths, usually using a master equation [301 EH |33l IMl EH EEl EZ] solving it by 
either numerical integration, diagonalization, or Monte Carlo wave-function approach. 
Using these methods one can study chains of up to n = 20 spins [3ll EZ] • 

Using MPO approach for solving master equation proposed in the present paper 
we can calculate NESS of ~ 100 spins. Even though we are at present unable to give a 
precise statement about the efficiency of the method for a general NESS (for an exception 
due to exact solvability see [E]), MPO approach should nevertheless prove very useful 
in studying out-of-equihbrium quantum many-body phenomena. 

In section |2] we outline our numerical algorithm of computation of NESS in terms 
of MPO ansatz as an asymptotic time-dependent solution of quantum (Lindblad) 
master equation. We also describe in details efficient models of the baths used in our 
calculations. In sections El and H] we discuss our numerical results for NESS and their 
transport properties in several integrable and non-integrable models of quantum spin 
chains, whereas in section [5] we summarize our findings and conclude. 

2. Master equation 

The most general completely positive trace preserving flow can be written in a form of 
the Lindblad equation [38] (see [39] for a comprehensive review on completely positive 
flows) 

= i[p, H]+^J2 (l^'^/^' 4] + [Lk, pLl]) , (1) 

k 

(setting ^ = 1) where H is the Hamiltonian, generating the unitary part of evolution, 
Lfc are Lindblad operators, and 7 is some overall bath-coupling strength parameter. 
Formally, the solution of linear Lindblad equation can be written as p{t) = exp (£t)p(0), 
where £ is a Liouvillean super-operator corresponding to the right-hand-side of the 
Lindblad eq.([T]). 

2.1. Matrix product operators and simulation of time- evolution 

We are going to simulate the time-evolution under the Lindblad equation by using a 
matrix product operator (MPO) formulation. For a chain of n spins 1/2, the density 
matrix p can always be expanded over all possible products of local Pauli operators 
which form a basis of 4" dimensional Hilbert space of operators. 
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where we use short notation a- = cr^^ ■ ■ ■ a*", s = Si, . . . s„, and Sj G {0, 1, 2, 3}, with 
(7° = 1,(7^ = cr^, o"^ = 0"'', = a^. A lower index in Pauh operators, for instance I in af , 
will always denote a position / G {1, ... n} of the spin on which it operates. In the MPO 
ansatz, the expansion coefficients Cg are represented in terms of An D x D matrices A*% 
i = 1, . . . , n, as 

c^= tr(Ar---A-). (3) 

The propagator exp (£t) is implemented within the MPO formalism in small time 
steps of length r (in our simulations r = 0.05). We decompose Liouvillean as 
C = Ci + £2, with the condition that all the terms grouped within each C^, z/ = 1,2, 
mutually commute. For example, £1 contains terms with interactions between 2nd 
and 3rd spin, 4th and 5th spin, and so on, while £2 contains interactions between 
other pairs, whereas the corresponding one-body terms are distributed evenly between 
£1 and £2. For each small time step we use Trotter-Suzuki formula, exp (£r) = 
Jlj.exp (Q!fc£ir) exp (/9fc£2T) + C(t^), with an appropriate scheme having either p = 3 
or p = 4, depending on the required accuracy. MPO representation of a density matrix 
p is exact only if the matrix dimension D of a matrices at l-th site is equal to the 
number of nonzero Schmidt coefficients fij for a bi-partite splitting of a super-ket |p) to 
the first / sites, denoted by A, and the remaining n — I sites, denoted by B, 

j 

with orthogonal vectors and \^f)- The inner product is defined by {a\P) = 
tr (q;^/9)/2". Since the maximal number of nonzero Schmidt coefficients is usually 
exponentially large in n one truncates exact representation by keeping only a small 
fraction of largest Hj. An error made in such a truncation is given by a sum of truncated 
[IQ]. As a rough estimate of a minimal needed matrix dimension D one can use Von 
Neumann entropy of a super-ket denoted by S^, also called operator space entanglement 
entropy (OSEE) (for an earlier definition in a different context see [41j). In terms of 
coefficients Cg OSEE can be expressed as, 

= -tr^(Rlog2R), R={p\p)-Htb\p){p\, (5) 

where, again, a subscript A denotes the first / sites and B its complement. Writing the 
expansion coefficients 2" dimensional matrix C^^^g^ = c^^s^, the matrix R 

(IS]) is given by R = CC"^/ tr (CC"^). OSEE is a measure of entanglement of a super-ket 
\p) which, however, is essentially different quantity as entanglement of a mixed state 
represented by the density operator p. Also, OSEE is different than the ordinary von 
Neumann entropy of p. For instance, OSEE either saturates or grows logarithmically 
with time in integrable transverse Ising model [12] or in Heisenberg XXZ model in a 
random magnetic field [13], it saturates or grows linearly with n in NESS of XY spin 
chain ^3], and it saturates or grows logarithmically with the inverse temperature for 
generic thermal states [18]. 
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2.2. Physics and implementation of the baths 

Let us say a few words about the details of our implementation of Markovian baths, 
that is of a non-unitary part Cb of C, involving Lindblad operators Lk. In our case - as 
our goal is to simulate many-body coherent transport in one dimension - the baths, i.e., 
the Lindblad operators L^, will act only at the two ends of a ID chain. More precisely, 
they will act just on the first and the last spin for a single-spin bath, or just on the first 
two and the last two spins for a two-spin bath. For a single spin bath, we put Cb as part 
of Ci as it commutes with all the other terms, so in the corresponding Suzuki- Trotter 
propagator exp (Cbt) factorizes out. Our choice of Lj. will be such that Cb will have a 
single non-degenerate eigenvalue equal to with the corresponding eigenvector pB being 
a local equilibrium state, i.e. a direct product of thermal states of (each of) the edge 
spins, or the pairs of the edge spins. Furthermore, due to stability of completely positive 
dynamics, all other eigenvalues of Cb should have negative real parts|il 

As a consequence, p{t) will for sufficiently long time converge to NESS poo = 
limi^oo exp (£t)p(0) of the entire spin chain, if starting from a generic initial state 
p(0) having a non- zero overlap with pac- Due to ergodicity of the internal coherent 
dynamics generated by the Hamiltonian if, we assume: (i) the resulting NESS poo will 
be unique, i.e. independent of the details of the initial condition p(0), (ii) the resulting 
NESS Poo will be (for complex - say non-integrable and quantum chaotic - internal 
dynamics) independent of the details of the bath operators - and will only depend 
on thermodynamic bath parameters, such as temperature, magnetization, chemical 
potential, etc. Both assumptions have been carefully checked and confirmed in the 
numerical simulations which are reported bellow. 

Since non-unitary evolution does not preserve Schmidt decomposition structure of 
matrices A^' we apply local rotations every few steps in order to "reorthogonalize" 
matrices A^', recovering the Schmidt decomposition [2]. We note that Ref. j9] originally 
proposed the idea to use MPO for time-dependent solutions of Lindblad master equation, 
however it implemented it in a physically essentially different context of systems with 
bulk-dissipation where each spin has been monitored by a Lindbladian bath. The fact 
that operator space entanglement is rather small in such a situation is not surprising, 
neither is the fact that it cannot describe coherent out-of-equilibrium many-body 
phenomena. 

The goal of the present paper is to test the efficiency of MPO simulation of NESS 
as given by the asymptotic t —>■ oo solution of the master equation driven only by 
the boundary Lindblad terms, for different (integrable and non-integrable) ID spin 
chains, checking along the way their transport behavior being for instance that of 
a normal (diffusive) conductor or displaying anomalous transport. One expects that 
quantum chaotic systems will display normal conduction [221 132] , that is, they will obey 
Fourier/Fick/Ohm's law in which the current j is proportional to the gradient of a 

if For an exact analysis of such open out-of-equilibrium quantum dynamics for Wigner-Jordan solvable 
models see Refs. [TSlfM] , 
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driving field e, say local energy density /temperature (or spin density /magnetization, or 
particle density/chemical potential, etc) 

j — —kVs. (6) 

For such normal conductors the transport properties should not depend on the details 
of the baths (assumption (ii) above). On the other hand, for the integrable models the 
choice of the baths could play a role. Wc arc going to use two different models of the 
baths. In studies of spin conduction in the Heisenberg model we will use a single-spin 
bath while we are going to use a two-spin bath when studying energy transport in a 
quantum chaotic tilted Ising model. The reason for using a two-spin bath was in a better 
(faster) convergence to NESS with time t and size n, which is probably due to the fact 
that the energy-density is a two-body operator, while the spin-density is a one-body 
operator. 



2.3. Single-spin hath 



When studying the spin (magnetic) transport we are going to couple the first and the 
last spin of the chain to a single-spin bath. For simplicity we shall bellow write only 
one part of Cb pertaining to a given edge spin, being either I — 1 or I — n and omitting 
the index I. In all our numerical simulations with the single-spin bath we shall take the 
coupling strength 7 = 1. 

There will be two Lindblad operators acting on the spin. 



1 



r± = 



r+(a^ + iay), 

' l =F tanh^L,R 
1 ± tanh/iL,R' 



r_((7"-i(7y) 



(7) 



Stationary state 

p-^diag(r+,r_; 



for £b constructed from the above operators is pe = 
oc exp (— /iL,RC"^), with the average magnetization tr (pB^"^) = 



— tanh /iL,R- We stress that this is a stationary state of a single spin in the absence 
of Hamiltonian evolution. Parameters pl,r of the left and right bath, respectively, play 
a role of an external thermodynamic potential enforcing a spin-density gradient, say a 
magnetization of a macroscopic magnet in contact with the edge spin. Matrix represen- 
tation of a super-propagator exp {C-qt) in the Pauli basis cr", a = 0, 1, 2, 3, reads 
/I \ 



exp {Cbt) 









-(r++r_)r 








-(r++r_)T 








g-2(r++r_)r y 



(8) 



2.4. Two-spin hath 

Here we would like to construct a bath i-e. determine the corresponding Lindblad 
operators L^, which produce a given unique stationary state pB of a pair of spins, 
e.g. such that pb = exp(— /i/T)/ tr exp(— /i/T) is a local thermal (Gibbs) state with 
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respect to some two-spin energy density operator h = /i((tJ^, cxg^). In all our numerical 
simulations with the two-spin bath we shall take the coupling strength 7 = 2. 

For the description of the methodology we shall assume that pe is a known but 
completely general state given in terms of a 4 x 4 matrix. We require that pb is a unique 
eigenvector of £b with the corresponding eigenvalue 0, while all the other eigenvalues 
of Cb are negative. We are looking for a set of such that the resulting £3 will have 
the above properties. With the constraints so far the choice of L^'s is not unique as we 
have fixed only one eigenvector and the corresponding eigenvalue of Cb- So in addition 
we shall require that all other (negative) eigenvalues of £3 are equal to —1. Such Cb 
will produce the fastest possible convergence to pe for a given fixed spectral norm of 

Assuming first that pe = diag((io, (ii, ^2, (is) is diagonal one can easily check that 
the following set of Lindblad operators Lk, 

Lij = \l^r'^r\ 2,j = 0,1,2,3, (9) 




m = {i mod 2) + 2 (j mod 2),r* = {ct^ + i a^, - i a^, 1 + ^7^ 1 - a"}, 

results in £3"^^ satisfying the above conditions. Note that the above 16 Lindblad 
operators, labeled by a double index k = (ij), can in fact be replaced by a set of 
15 traceless operators leading to the same £3"^^. In the Pauli basis the only nonzero 
matrix elements (^b'^^)^,/?, a, j3 E {0, . . . , 15} are (assuming positive di and tr pe = 1) 

{^t^^)a,a = - 1, a = 1, . . . , 15 
(£g'^^)i5,o = do — di — d2 + (is, 
(-^B )i2,o = do + di — d2 — (is, 

(£f^)s,o =do-d^ + d2-ds. (10) 

Basis states are enumerated in such a way that the least significant bit is the first one, i.e. 
corresponding to the left factor in tensor products Qj. Then, to obtain the bath data 
Lk and Cb for a non-diagonal pb we write the eigenvalue decomposition of pb = V^d V, 
where d is diagonal and V is unitary. In terms of the diagonal part d we first obtain 
the Lindblad super-operator £3"^^ as described above ( iQllTOl) and then rotate it in the 
operator-space using the transformation R induced by V. Writing the orthogonal matrix 
of R in the Pauli basis a- = o""^ ®cr"2 -we have -Rq,/? = tr (V^V- V"(T^)/4, giving the final 
Lindblad propagator 

exp (£b r) = R^ exp {C'^''^ t)R. (11) 



3. Spin transport 

Here we are going to study the spin transport in Heisenberg XXZ model with the 
Hamiltonian 

n— 1 n 

(=1 1=1 
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Figure 1. Spin profile (a) and spin current (b) for Heisenberg XXZ model at 
A = 0.5 and bath couplings (0 with /xl.r = ±0.22. Data for different chain lengths 
n = 16,20,32,64, from green (bright) to black curve, are shown. With increasing n 
the spin profile gets flat - an indication of an ideal spin conduction, visible also in the 
independence of the average spin current j on the chain length n. 



The first and the last spin will be coupled to a single-spin bath (JHl). The initial state is 
chosen to be a product state p(t = 0) oc exp (— Yliil^i^l)^ where /i; linearly interpolates 
between the left/right bath values /xl,r- We find that for times t of the order of several 
times n the state practically converges to NESS poo, the properties of which we are 
interested in. In particular, in order to asses the validity of the spin Pick's law we 
are going to calculate the magnetization profile (also referred to as the spin profile), 
Si = tr {afp^), and the local spin current defined as 

ji = tr [(o-r<i - ^r+i)poo] . (13) 

We are going to study three different parameter regimes of the Heisenberg model 
corresponding to qualitatively different nature of many-body dynamics. The integrable 
XXZ Heisenberg model in the absence of magnetic field is known to display an ideal 
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Figure 2. Spectrum of Schmidt coefficients /i^ for a symmetric bipartite splitting and 
different MPO dimensions D = 40, 80, f 20 of a NESS in the Heisenberg model ^ 
with A ~ 0.5 coupled to two single-spin baths with /il,r = ±0.22. Full (red) curves 
are for the chain size n = 64 while dashed (blue) curves are for n = 20. Dotted line 
indicates the exponential decay exp(— 1/40). Inset shows the dependence of OSEE 
S"" ([5]) on n. does not grow appreciably with n, indicating the efficiency of MPO 
representation of NESS. 

spin conduction for A < 1 while it is probably a diffusive (normal) spin conductor 
for A > 1 pni EH EHl Elj. The last statement is quite controversial in the light of 
the algebraic integrability of the model [H] so it will be inspected more carefully in 
the present paper. As the third case we shall study spin transport in the XXZ model 
with staggered transverse magnetic field in the regime of quantum chaos, significantly 
improving numerical evidence for the conjecture (put forward for the case of heat 
transport in [221 [22]) that quantum chaos corresponds to normal diffusive transport 
(Fourier's, Ohm's or Pick's law). 

3.1. Ideally conducting regime, A = 0.5 

We will first take A = 0.5 and no magnetic field, hi = 0, in order to test the method in a 
regime where Heisenberg model is an ideal spin conductor (having a non- vanishing Drude 
weight at any, say infinite temperature [21]). The left and the right bath parameters 
are set to /iL = 0.22 and /iR = —0.22. The results of numerical simulation are shown 
in figure [TJ One can clearly see that the system is an ideal conductor, that is, it 
does not obey spin Pick's law because the spin current is found proportional to the 
magnetization difference and not its gradient. This is a typical property of completely 
integrable systems. 
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Considering numerical efficiency of MPO simulation of NESS, we show in figure [2] 
the Schmidt spectrum for a symmetric bipartite cut of NESS and different MPO 
dimensions D. Numerical data suggest that the tails of the spectra of Schmidt 
coefficients decay exponentially. OSEE ([5]) does not seem to grow with n (although 
based on numerics we can not exclude a very slow growth) suggesting the method is 
efficient, i.e. computation time is asymptotically polynomial in n (linear in n if OSEE 
converges) . 

3.2. Normal conductor, A = 1.5 

Increasing A above 1 the Heisenberg model becomes a normal (diffusive) spin conductor 
(but not a diffusive heat conductor!) despite its integrability. We take A = 1.5 in the 
absence of the magnetic field, hi = 0, and a single-spin bath with a weak driving 
/iL,R = ±0.02 ([8]) in order to make sure that non-linear transport features do not 
obscure the effect. Figure [3] clearly demonstrates that the system is indeed normal 
spin conductor. Spin profiles in NESS are, apart from boundary effects, perfectly 
linear. To reduce the boundary effects we drop the leftmost and the rightmost two 
spins when calculating the drop of magnetization AS* = Sn-2 — S3, and its gradient 
VS = AS/{n — 4:). Spin current is clearly proportional to the gradient of magnetization 
VS, or at fixed bath data, to ~ 1/n, see figure Eb. Note that the largest chain size 
n = 100 is much larger than what has been numerically achievable with other methods, 
like various Monte Carlo wave-function methods [32l [Ml 135] . 

To show the convergence of p{t) to the asymptotic NESS poo we show in figure H] 
the snapshots at different instants of time t of local spin current proffies ji ( |T3l) . Clearly, 
a tendency towards a uniform current profile at NESS required by continuity equation 
is observed. 

The tail of the Schmidt spectrum decays here qualitatively slower than for an 
ideally conducting case (e.g. A = 0.5 of fig. Wp), exhibiting perhaps an asymptotic 
power law decay /ij ~ l/i^-"^^. Note that if Schmidt coefficients decay algebraically 
as /ij ~ l/i^ with some power p, then the OSEE S*" converges to its exact value for 
D ^ 00 as Sjj^^ — S*!) ~ log2 (D) / D"^^'^ . For our p ^ 1.25 this would mean a very slow 
~ log2 {D)/D^'^ convergence of SK Since we are limited to relatively small dimensions, 
of order D ~ 100, we are here not able to asses whether 5*^ saturates with n or not. 

3.3. Normal conductor, A = 0.5 and staggered transverse field 

Here we are going to take the Heisenberg model (fT2l) with A = 0.5 and staggered 
transverse magnetic field, h2i+i = and /i2/ = —1/2. We use the bath parameters 
/iL,R = ±0.1. We have checked that for these parameter values the spectrum of the 
Hamiltonian f|T2l) exhibits the characteristics of quantum chaos, i.e. the energy level 
spacing distribution agrees with universal predictions of random matrix theory with 
no free parameters, so the spin conduction is expected to be normally diffusive and 
to obey the Fick's law. The results shown in Figl6] clearly conffim this expectation, 
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Figure 3. Spin profiles and current versus system size in NESS for the Heisenberg XXZ 
model at A = 1.5 and the bath coupling parameters /iL.R = ±0.02. Data for different 
chain lengths n = 16, 32, 64, 100, from green (bright) to black curve, almost overlap 
(frame (a)). Apart from the boundary effects the spin profile is linear, suggesting 
normal conduction and diffusive transport. This is confirmed in frame (b), where we 
show the dependence of the scaled spin current j /AS on the chain length n. Dotted 
line is j / AS 1.15/(n — 4), indicating spin "Pick's law" with the spin conductivity 
K = 1.15. 



namely the spin density profiles are linear, and the spin current decays as oc 1/n. Note 
that this figure is very similar to Figj3]for the diffusive-integrable case with A = 1.5, a 
subtle difference may be that now the boundary effects seem to be stronger. We have 
thus dropped border 5 spins at each end of a chain when calculating the spin drop 
AS = Sn-5 — Sq and the gradient VS = AS/ {n — 10). 

The tail of the Schmidt spectrum shown in figure [7] seems to decay even slower now, 
/ij ~ l/i'^'®°, making the simulation harder than in both integrable cases without the 
magnetic field. This is consistent with the previous observations in MPO simulations of 
long-time evolution in the Heisenberg picture [11] , 
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Figure 4. Convergence of the spin current profiles (dependence of the local current ji 
on the lattice site i) with time for the Heisenberg XXZ model at A = f .5 (the case of 
Fig. [3]) and for size n = 40. Snapshots at times t = 4, 10, 20, 40, 80, 170 are shown, for, 
respectively, automatically adapted MPO dimensions D — 40, 60 and later D = 80. 
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Figure 5. Spectrum of Schmidt coefficients /i^ of NESS (at convergence time t — 250) 
for a symmetric bipartite cut of the Heisenberg XXZ model at A = 1.5. MPO 
Dimension is D — 120 and the spin chain length n ~ 64. Dotted line indicates a 
power law decay with finite size effects setting in when i ^ D. 
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Figure 6. Spin profile (frame (a)) for quantum cliaotic Heisenberg XXZ model at 
A = 0.5 in a staggered transverse field hi — (0,-0.5,0,-0.5,...). Data for chain 
lengths n — 20, 30, 40, 80, from green (bright) to black curve, are shown. Apart 
from the boundary effects the spin profile is linear, suggesting normal diffusive spin 
conduction. This is confirmed in frame (b), where we show the dependence of the 
scaled average spin current j'/AS* on n. Dotted line is j/AS ^ f 8.0/ (n — f 0), indicating 
"Pick's law" with the spin conductivity k — 18.0. 



4. Energy transport 

As the last test we are going to study the heat conductivity in Ising model placed in a 
tilted magnetic field, tilted Ising modelioi short, with the Hamiltonian 

n-l ^ ^ 

H = Y^ /iM+i = -2fTfaf+i + -{Kaf + + -(/i^+i + ^z^xf+i), (14) 

1=1 

with /ix = 3.375 and /i^ = 2. For these parameters the system is quantum chaotic, see 
e.g. [32j|, and displays normal heat conduction [32l[3l], however, previous simulations 
were limited to system sizes n < 20. For studies of heat transport in other quantum 
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Figure 7. Spectrum of Schmidt coefficients of NESS (at convergence time t = 200) 
for a symmetric bipartite cut of quantum chaotic Heisenberg model in a staggered 
field. MPO Dimension is _D = 150 and size n = 50. Dotted line indicates a power law 
decay 



chaotic systems see also [221 ESj. Local energy current is in this case 

ji = 2K tr [i(Tf_,aJ - afaf^,)p^] . (15) 

When we used the single-spin bath (|H]) the convergence of local energy current to a 
homogeneous site-independent value expected for NESS was rather slow, meaning that 
the spectral gap of the quantum Liouville operator C is inconveniently small for such a 
bath model. Therefore, we rather used a two-spin bath ( ITTi) for which these effects are 
smaller. We set the temperature of the left bath to Tl = 20 and of the right bath to 

Tr = 30. 

Since in an out-of-equilibrium system the definition of local temperature may not 
be completely unambiguous we have instead of looking at the temperature profile, 
computed the energy density profile ei = trpoo^«,«+i- As discussed in [32], the energy 
density uniquely determines the temperature in the equilibrium (thermal) state and it 
can also be used as a good measure of local temperature out of (but not too far from) 
equilibrium. 

The boundary effects seen in figure El are again relatively strong. Because of these, 
and because we kept Tl^r - which is a "non-interacting temperature" of a 2 spin system 
- fixed for all n, the actual local temperature varies slightly with n. In figure [Hk we 
therefore appropriately shifted individual energy profiles in vertical direction in order to 
get scaled overlapping curves. To decrease finite-size effects at boundaries we dropped 
border 4 spins when calculating the energy difference AE or the energy density gradient 
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VE = AE/{n-9) (fig.[H]). 

Taking all these into account, we have confirmed the Fourier's law of heat 
conduction j = K,AE/{n — 9), where AE was the energy density difference between 
spins (5, 6) and (n — 5, n — 4), with an excellent numerical accuracy for 10 < n < 100. 

We have also looked at the tail of the spectrum of the Schmidt coefficients which 
again exhibits power law decay, best fitted with /Xj oc 




Figure 8. Local energy density profiles in NESS for quantum chaotic tilted Ising 
model ((HI). Data for chain lengths n = 20, 32, 50, 64, 100, from green (bright) to black 
curve, almost overlap (frame (a)). Apart from the boundary effects the energy density 
profile is linear, suggesting normal heat conduction. This is confirmed inframe (b), 
where we show the dependence of the scaled energy current j /AE on n. Dotted line is 
j / A.E ^ 5.5/(n — 9), indicating the Fourier law behavior with the conductivity k — 5.5. 



5. Summary and discussion 

In the present paper we have demonstrated numerical simulation of NESS for strongly 
but locally interacting open quantum systems in ID in terms of the MPO ansatz. We 
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Figure 9. Spectrum of Schmidt coefficients /i,; of NESS (at convergence time t = 500) 
for a symmetric bipartite cut of quantum chaotic titled Ising model (for the case of 
fig. [8]). MPO dimension is D = 150 and size n — 50. Dotted line indicates a power law 
decay i^^-'^'^. 

have considered the most difficult case where dynamics in the bulk is fully coherent 
(Hamiltonian), and dissipation (coupling to the baths) only takes place at the boundary 
(ends of the chain). In this setting we have been able to con&m the laws of diffusive 
transport, such as Fourier's law of heat conduction and (spin) Pick's law for spin 
conduction, in the cases where the underlying model is strongly non-integrable and 
displays the features of quantum chaos, or even when it is integrable but all the conserved 
quantities are irrelevant for the transporting current, like in the case of spin-conduction 
in Ising-like Heisenberg XXZ chain. The main purpose of our paper was to demonstrate 
that such simulations were now possible for system sizes of order of n = 100 spins 
1/2, which is considerably larger than with the competing methods (usually based on 
Monte Carlo wave-function techniques, where presently only n ~ 20 is reachable for 
fully out-of-equilibrium simulations). 

As for the quantitative analysis of the efficiency of MPO ansatz for NESS we 
have analyzed the spectrum of the Schmidt decomposition for the worst-case (half-half) 
bipartition of the chain, when treating the density matrix of NESS as an element of the 
Hilbert space of operators. Summarizing these results, we have found that completely 
integrable and ideally conducting cases, in our example XF-like Heisenberg XXZ chain 
(|A| < 1), are easiest to simulate since there the tail of the spectrum of Schmidt 
coefficients decays exponentially so MPO of a rather modest matrix dimension gives 
a good description of NESS. 
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For still integrable systems, but for which all the conservation laws are irrelevant to 
the transporting current in NESS, such as the Ising-like Heisenberg XXZ chain (| A| > 1), 
we have found qualitatively worse efficiency, namely there the tails of the Schmidt 
spectrum appear to exhibit power law tails fii cx i~^-'^^. 

At last, for non- integrable systems in the regime of quantum-chaos, the efficiency of 
simulation appears again worse than in the integrable cases above (in our examples we 
have looked into the Heisenberg XXZ chain in a staggered transverse field, and Ising spin 
chain in a tilted magnetic field). Namely, there we have found algebraically decaying 
tails of the Schmidt spectrum /ij oc with the power p G [0.7, 0.8]. 

Based on our results we also conclude that zero-temperature (ground state) 
properties of the system, whether being critical or gaped, do not infiuence far-from- 
equlibrium properties of NESS. 

The work has been financially supported by Program Pl-0044 and Grant Jl-7347 
of Slovenian Research Agency (ARRS). 
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